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ABSTRACT 

We discuss chemical reaction networks and metabolic pathways 
based on stoichiometric network analysis, and introduce deformed 
toric ideal constraints by the algebraic geometrical approach. With 
the deformed toric ideal constraints, the shape of flux is constrained 
without introducing ad hoc constraints. To illustrate the effective- 
ness of such constraints, we discuss two examples of chemical re- 
action network and metabolic pathway; in the former the shape of 
flux is constrained completely by deformed toric ideal constraints, 
and in the latter, it is shown the deformed toric ideal constrains the 
parameters of flux at least partially. 
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1 Introduction 



The metabolic pathway analysis is one of the major fields in the systems 
biology and is the basis of the metabolic engineering and synthetic biology. 
With the development of experimental technology, this field is studied with 
practical objectives, such as clarifying the metabolic systems of life with bio- 
logical interests and manufacturing biochemical products through metabolic 
processes. 

In chemistry, it is important to identify the chemical mechanism. The 
instability of steady states of chemical reaction systems exhibits exotic dy- 
namics, such as switching between multiple steady states, explosions, and 
sustained oscillations. The chemical reaction networks have been studied 
with the purpose of clarifying the chemical mechanism behind such phenom- 
ena. 

Theoretical studies on chemical reaction networks are based on stoichio- 
metric network analysis (SNA), which is the method based on the mass action 
kinetics. SNA was initiated by Clarke [H |2l |3] and succeeded by Feinberg 
et al. [H El El [7] . Gatermann et al. [HI El [10] have studied chemical reaction 
networks from the viewpoint of algebraic geometry, especially of polynomial 
rings. Inspired by their study, Shiu et al. have used toric varieties in the 
analysis [HI |12l |13]. There are also studies on bifurcations of dynamical 
systems with SNA [H [ig [161 [IT] . 

Metabolic pathways, on the other hand, have been studied by fiux balance 
analysis (FBA) which introduces the steady state solution space, called 'null 
space' or 'kernel space.' The introduction of 'elementary modes' [T8 l \T9 \ [20 l 
[2TI [22] and 'extreme pathway' [23] is specific progresses in this field. Review 
articles on these theme include refs. [2ll[25||. For the discussions of these two 
approaches, see ref. [26] . 

This paper considers a unified approach to obtain constraints on the fiux 
for both fields of chemical reaction networks and metabolic pathways, and 
applies SNA to analysis of metabolic pathways. While the theoretical anal- 
ysis in these two fields have common bases, the steady state solutions have 
been analyzed in different ways, and there have been few works discussing 
them in a unified way. Ref. [23] is an exception, where the relationships be- 
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tween SNA and FBA are discussed, but they do not apply SNA to analyze 
the metabolic pathways. We discuss algebraic geometrical constraints on sto- 
ichiometric networks: from the monomial vector expression of flux, we derive 
the deformed toric ideal, that works as constraints on the parameters of the 
flux. The existence of such constraints had not been pointed out before. 

We show two examples of analysis by the proposed method in a chemical 
reaction network and a metabolic pathway. In the former example, we show 
that all of flux parameters are determined completely by the reaction coeffi- 
cients with the deformed toric ideal constraints. In the latter, we show that 
the flux parameters are partially constrained by the deformed toric ideal, 
without which the flux parameters cannot be restricted. We also describe 
the relation between the deformed toric ideal constraints and thermodynamic 
feasibility. 

This paper is organized as follows. In section 2, we review the basic 
formulation of stoichiometric network analysis with an example of chemical 
reaction network. In section 3, we explain stoichiometric flux. In section 4, 
we discuss deformed toric ideal constraints. In section 5, we apply the argu- 
ments of SNA and deformed toric ideal contraint to an example of metabolic 
pathways. In section 6, we give the conclusions of the analysis of the current 
paper and the future directions of mathematical studies of stoichiometric 
networks. 

2 Stoichiometric Network Analysis 

The stoichiometric network analysis (SNA) starts with the chemical reaction 
systems which are described by 



Here, x and x are a concentration vector of reactant species and its time 
derivative, respectively, N is the stoichiometric matrix, and v{x\ k) is the flux 
vector, where k is reaction coefficients. The stoichiometric matrix, explained 
below, can be determined, once we know the form of chemical equations. 



X — N ■ v{x] k). 



(1) 



aijSi -|- • • • -|- djnjS, 
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bijSi + ■ h bmjSm, j — ^, ■ ■ ■ ,1, 
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where kj is the reaction coefficient for reaction j. The form of flux vector 
v{x] k) is determined by the mass action kinetics. By the mass action kinetics, 
the verocity of reaction is described by 

kASir^-'-lSmT-^ (3) 

where [Si] denotes the concentration of Si. With the matrices A and B of 
the elements aij and bij, respectively, is deflned as 

N = B-A (4) 

The stoichiometric equation ([1]) shows that the time derivative of concen- 
trations of reactant species is represented by the product of the stoichiometric 
matrix and flux vector. 

Example 2.1(Gatermann et. al. 2005 |[14|) 

We consider the phosphofructokinase reaction which is a part of glycosis. 
It is an extension of a reaction system proposed by Sel'kov '68 . There are 
m = 3 chemical species; 5*1 denotes the product Fructose- 1,6-biphosphate, 5*2 
denotes the reactant Fructose-6-phosphate and the extension 5*3 stands for 
another intermediate in equilibrium with Fructose-6-phosphate. The / = 7 
reaction laws are given by 

2Si + S2^ 3Si 

kd ks kj 

52^o45i453. 

We arrange the left and right hand side of reaction laws ([2]) in so-called 
complexes Cj, j = 1, . . . ,n with the rate constants kij > 0; 

^56 ^64 ^43 

L'5 ^ L'6 ^ L'4 ^ tv3- 

kes kie ^34 

By the mass action kinetics, this reaction system is described by the 
following differential equations. 



Xi — k2ix\x2 + /C46 ~ k^/^Xi — /i;34Xi + ^43X3, 
X2 = -k2ix\x2 + k^(^ - kQ^X2, 
= k^^Xi — k43X3. 
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Represented by the stoichiometric equation form, and v{x; k) in eq.(IT]) are 
given by 



N 



11-10 0-1 
-10 1-10 
1 



1 





1 



) 



(5) 



and 



/ k2ixjx2 \ 



v(x; k) 



km 
k^Xi 

^56 
k34Xi 



(6) 



V ^433:3 / 



3 Stoichiometric Flux 



The steady state solution for a chemical reaction (dynamical) system is de- 
fined by the concentrations whose time derivatives vanish (a; = 0), i.e. 



where z = v{x]k). When there is no confusion, we use z hereafter. The 
solution of z that satisfies eq.© forms a convex polyhedral cone, and the 
steady state solution space is called 'kernel space' or 'null space.' Minimal 
generating vectors of the steady state solution are called 'extreme currents.' 
Each extreme current is the generator of a convex polyhedral cone. A non- 
negative linear combination of the extreme currents is also a steady state 
solution. 

Example 2.1 (continued) 

The extreme currents satisfying eq.([7]) are given by 



Nv{x; k)^Nz = 0, 



(7) 
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If we take the linear combination of the extreme currents, 



JlEl+j2^2+j3^3+j4^4 
/ J4 \ 

3i 

Jl + J4 
J2 + J4 

h 

J3 



v 



(9) 



(10) 



where, ji are nonnegative hnear combination coefficients. Thus, 2; is a general 
steady state solution of eq.Q. While the coefficient ji might look to take an 
arbitrary nonnegative value, they have some constraints as will be shown in 
the next section. 



4 Deformed Toric Ideal Constraints on Stoi- 
chiometric Network 

By the explicit monomial vector form, the steady state solution forms a 
deformed toric ideal. This can be seen by describing the elements of z as 
monomials in the original x coordinates (in affine space). The ideal ly/ = 
{f & Q[z]\f{x) =0} C Q[2;] is called a deformed toric ideal, where Yl is the 
configuration whose entries are the exponents of the monomials in the flux 
vector. 
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We will see that the form of flux is constrained by the deformed toric 
ideal. For the mathematical basics, see refs. |28l |29] and for the derivation 
of toric ideals, see ref. [30] . 
Example 2.1 (continued) 

For the current example, the generators of the deformed toric ideal are 
obtained as a binomial form, and it is easy to calculate from eq.([n]); 

Iz = {^2 — ^46) ^21^5^3 ~ ^64^65^1) ^34^3 ~ ^64^6) ^4 ~ ^Se)- (H) 

By the correspondence between the elements of vector z in the repre- 
sentation by zi and by j/, the deformed toric ideal can be described by ji 
coordinates. By substituting ffTOj) with f|TT]) . we obtain deformed toric ideal 
in ji coordinates. 

= O'l - ^46, A;2ij2(ji +3if - klih^ii, hih - +ii),i2 +3i - he) (12) 

In the rest of this section, we show that ji are constrained by the deformed 
toric ideal. As a result, ji are determined without introducing ad hoc con- 
straints. Note that in the following derivation we use only the generators of 
deformed toric ideal, that is, the relationships between monomials. This has 
not been used in any previous studies on this example. From the constraints 
that each generators of deformed toric ideal are zero, we obtain the following 
expressions of ji , 

ji = he, 

ja = -Z + ^56, 

is = k3i{k^Q + Z) / kei, (13) 
34 = Z. 

Here, Z is the solution of the following 3-dimensional algebraic equation, 

k2lZ^ + {~k2lk56 + 2k2lki6)Z^ 

+ ik2ikle + klJtQ^ - 2^21^46^56)^ - hikjek^Q = 0. (14) 
The positivity conditions for ji give 

0<Z < k56. (15) 
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Eq.f lTS]) shows that the hnear coefficients ji depend on the value of reaction 
rates kij and Z. Eq.f lT4|) shows that Z also depends on the value of kij, 
and thus ji depends on the value of kij, which means ji cannot be chosen 
arbitrarily, indifferent to the value of k^j. Note that these constraints are 
derived by the algebraic property of monomials, without introducing any ad 
hoc constraints. This strong constraint seen in the current example does not 
hold for general reactions: the example has many generators in the deformed 
toric ideal enough to constrain all ji by k^j. In general cases, ji are only 
partially constrained by k^j. 

Substituting ji, . . . , ^4 of eqs. (lT3|) for eq. (fTO|) . we obtain the stoichiometric 
flux z under deformed toric ideal constraints. 



(16) 



Z 

km 
km + Z 

^56 

-Z + k^Q 
kM{km + Z)/kQi 
\ ku{kA& + Z)/kQi J 

This means that we can determine the form of flux (length of extreme cur- 
rents), once we know the value of reaction coefficients kij by means of exper- 
imental observation. 



5 Example of Metabolic pathway 

In the above sections, we considered a chemical reaction network as an 
example. The same argument holds for metabolic pathways. In this sec- 
tion, we consider deformed toric ideal constraints with a concrete example of 
metabolic pathway. 

Example 5.1: 

Feedback Inhibition of pathway, Palsson (2011) |31] 

We use this example because it is one of the simplest realistic pathways whose 
monomial vector form of the flux is explicitly known. 
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Figure 1 : Feedback Inhibition of pathway 

In a biosynthetic pathway, the first reaction is often inhibited by the end 
product of the pathway. Fig{T] illustrates a protypical feedback loop in a 
biosynthetic pathway. 

A metabolic intermediate xi is formed and degraded as 



;i7) 



Then, if an enzyme Xq is expressed, Xi can be converted to X2'. 



Xi + Xq ^ X2 + Xq 

which is followed by a series of reactions 



(18) 



fc2 fca k4 fcs 

3^2 — ^ 2:3 — )■ X4 — )■ X5 — )■ 



(19) 



to form X5, the end product of the pathway. The end product has inhibitory 
feedback to the enzyme xq by binding to it and converting it into an inactive 
form: 



X6+ X^ ^ X7 
k-6 



(20) 



This system represents a simple negative feedback loop. The differential 
equations that describe this feedback loop are 
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Xi = 61 — koXi — kiXQXi, (21) 

X2 = hxeXi - ^2^2, (22) 

X3 = k2X2 - k3X3, (23) 

0:4 = k^xs — k^Xi, (24) 

x^ = k^Xi - k^x^ - {kex^XQ - k^QXj), (25) 

Xq = -kex^xe + k.^xj, (26) 

x^ = k^x^XQ — k^QX-j. (27) 



In the above equations, RHS is the sum of reaction rates. In ref . [31] (chapter 
2.2), reaction rates are described mathematically using kinetic theory. He 
discusses mass action kinetics as one of the fundamental concept of kinetic 
theory. 

We analyze this pathway with SNA, which was not done in ref. [3T]. In 
the stoichiometric matrix form, these differential equations are written as 



X = Nz, (28) 
and the metabolic pathway analysis starts with the following equation, 

Nz = 0. (29) 

The LHS is the product of the stoichiometric matrix and the flux vector, and 
we will find the steady state solution. This is the flux balance equation. For 
flux balance analysis (FBA), see refs. [32l[33], for example. 
For the current example, the stoichiometric matrix is 
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(30) 



and the flux vector is 
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/ kiXQXi \ 
koXi 

k2X2 

k^xs . (31) 
k^Xi 
k^x^ 

k-6X7 

\ h J 

For the metabolic pathway analysis, extreme currents correspond to ex- 
treme pathways [23]. The extreme pathways computed from the stoichio- 
metric matrix are 



El 



/0\ 
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VI/ 



( 1 \ 
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V / 



E^ 



/0\ 

1 






1 



(32) 



By taking a linear combination of the extreme pathways with nonnegative 
coefficients, the metabolic flux is obtained; 



3\E\ + J2^2 + J3^3 
/ 32 \ 
jl + J3 

3i - 32 

32 
32 
32 
32 

3i + is 
V 31 ) 



(33) 
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Here, we derive the deformed toric ideal of this pathway. Admitting Laurent 
monomials, from the monomial vector representation of z, the deformed toric 
ideal is given by 



h = {zikeko/zs - Z2kik^/zj). 



(34) 



From the corresponding representation of flux, the deformed toric ideal rep- 
resented by ii is given by 



There is only one deformed toric ideal constraint, which is obtained by 
equating the generator to zero. This is a relation between the extreme path- 
ways. The parameter region is partially constrained by the reaction coeffi- 
cients. In metabolic pathways, the shape of flux (length of extreme pathway) 
is important. The shape of flux is partially constrained by experimental ob- 
servations, furthermore automatically constrained by the deformed toric ideal 
constraints. 

In the above example, there is only one deformed toric ideal constraint. 
Comparing with the example of chemical reaction network, the number of 
constraints is small. This is caused by the small number of the generators of 
the deformed toric ideal. The smallness originates in the limited number of 
species which appear as the same form in the chemical equations. 

The points of ji = j^, 32 = or ji + ja = are critical points between the 
thermodynamically feasible and infeasible region in the sense of refs. [Mj 135], 
where thermodynamic feasibility is described as the orthogornality to the 
cycles of the oriented matroid associated with stoichiometric matrix A^. At 
the above points, the elements of flux is zero. For the detail, see ref. [36j- 




In the rest, we will show that, in addition to the steady state equations 
and mass balance, which were already known in ref. [31], the above deformed 
toric ideal constraint provides the complete solution to the steady state of 
this system. First, as in [31], the steady state equations are given by 



Ij = {j2k6ko/{ji -32) - j3)klk5/j2). 



(35) 
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= 6i — koXi — kiXQXi, (36) 

= kiXQXi — k2X2, (37) 

= k2X2 - ksx^, (38) 

= k^xs — k^Xi, (39) 

= ^4X4 — /C5X5 — {k^x^x^ — k^QXj), (40) 

= -kex^xe + k.^Xy, (41) 

= k^x^XQ — k_QXj. (42) 

It is known that these equations can be combined to give a quadratic equa- 
tion, 

y-^ + ay-b = 0, (43) 

where 



y = k2X2,a = k^k^Q^l + kiet/ko)/kQ,b = k^k^ekietbi/ (keko), (44) 

that has one positive root in y. Note that from y = ^2X2 = Z4 = j2 the 
variable j2 is given by a positive root of eq. pSI) . Next, the sum of Xq and x-j 
gives us the mass balance in the enzyme: 

xq + xj = et, (45) 

where et is the total amount of enzyme (constant). 

In ref . [31] , these constraints are not studied with ji coordinates. With the 
representation of flux by the extreme pathways, eq.( l45i) can be represented 
by ji coordinates, eq. (H5|) gives 

(ji + J3)h/j2k6 + (ji + j3)/k.e = et. (46) 

With these conditions, the region taken by the parameters ji is further lim- 
ited, combined with the deformed toric ideal constraint and nonnegativity 
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conditions on j/. 



ji = j2{hkQ + kih)/hh, (47) 
. _ 32{hk-Q{etkik(i - fcifcg - k^kp) - j2{kihh + fcjfeo)} ^ q 

Since from eq. psj) j2 is determined uniquely, eq.f H7|) shows ji is also unique. 
In addition, when the RHS of eq. fHHl) is nonnegative, a unique solution for 
(ji)j25 js) is determined, and when the RHS is negative, there is no steady 
state solution. 



6 Conclusions 

In this paper, we considered deformed toric ideal constraints on stoichiomet- 
ric networks, treating the chemical reaction networks and metabolic pathways 
in a unified way. This paper is the first that pointed out that the deformed 
toric ideal constrains the parameters of flux. We clarified that even the shape 
of flux is constrained by the deformed toric ideal. 

In general, for metabolic pathways, reactants and products are so diverse 
that chemical species may not appear as the same form in the LHS and 
RHS of the chemical equations. The number of generators for the deformed 
toric ideal is thus small and the effect of constraints may be limited. By 
considering sub-networks of large scale metabolic pathways, however, they 
should be more influential. It will be interesting to consider sub-networks 
like ATP production in energetic metabolism, because there appear some 
chemical species as the same form in the LHS and RHS of the chemical 
equations. 

In refs. [311 [35], thermodynamic constraints are considered explicitly, in 
terms of the non-equilibrium thermodynamic systems. Such constraints are 
not discussed in the current paper, but the study with thermodynamic con- 
straints, added to the deformed toric ideal constraints, will be interesting. 
By considering these constraints, flux is constrained in another form. 

As another topic, which is pointed out for the example of chemical reac- 
tion network, flux forms a convex polyhedral cone. The algebraic geometrical 
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study, using commutative algebra and combinatorics, will also be interesting. 
This is within our future work. 
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